function st=gen_s(T,n,xstar,ht,q,mu,v2)
st=zeros(n,T);
for tt=1:T
    for jj=1:n
        temp=q.*normpdf(xstar(jj,tt),...
            2*ht(jj,tt)+mu-1.2704,sqrt(v2));
        temp=cumsum(temp/sum(temp));
        st(jj,tt)=sum(temp<rand())+1;
    end
end
